analysis richness - repeated-measures of temporal data

library(nlme)
library(lme4)

## Loading required package: Matrix

## 
## Attaching package: 'lme4'

## The following object is masked from 'package:nlme':
## 
##     lmList

library(ggplot2)
library(sjPlot)
library(AER) 

## Loading required package: car

## Loading required package: carData

## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4

## Loading required package: lmtest

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: sandwich

## Loading required package: survival

library(brms)

## Loading required package: Rcpp

## Loading 'brms' package (version 2.12.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').

## 
## Attaching package: 'brms'

## The following object is masked from 'package:survival':
## 
##     kidney

## The following object is masked from 'package:lme4':
## 
##     ngrps

## The following object is masked from 'package:stats':
## 
##     ar

1. only for sea water treatments

1.1 load dataset (2 soil types x 4 inundation frequency by sea water)

# read dataset
df <- read.csv("alpha_diversity_dispersal.csv", row.names = 1, check.names = 0)
# only keep richness
df <- df[, -c(2:7)]

# remove samples from day 0
df <- subset(df, !grepl('___', rownames(df)))

# select soil from sea water treatment
df <- df[df$Water == "sea",]
df <- df[, -4]
df$Soil <- factor(df$Soil, levels = c("0", "70"))
df$Frequency <- factor(df$Frequency)
df$box <- factor(paste(df$Soil, df$Frequency, df$replicates, sep = "_"))
str(df)

## 'data.frame':    144 obs. of  6 variables:
##  $ richness  : int  828 1075 1101 1100 1108 1124 792 900 727 1074 ...
##  $ Soil      : Factor w/ 2 levels "0","70": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Frequency : Factor w/ 4 levels "f1","f2","f3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Day       : int  12 12 12 16 16 16 2 2 2 20 ...
##  $ replicates: Factor w/ 3 levels "A","B","C": 1 2 3 1 2 3 1 2 3 1 ...
##  $ box       : Factor w/ 24 levels "0_f1_A","0_f1_B",..: 1 2 3 1 2 3 1 2 3 1 ...

head(df)

##               richness Soil Frequency Day replicates    box
## 0_f1_sea_12_A      828    0        f1  12          A 0_f1_A
## 0_f1_sea_12_B     1075    0        f1  12          B 0_f1_B
## 0_f1_sea_12_C     1101    0        f1  12          C 0_f1_C
## 0_f1_sea_16_A     1100    0        f1  16          A 0_f1_A
## 0_f1_sea_16_B     1108    0        f1  16          B 0_f1_B
## 0_f1_sea_16_C     1124    0        f1  16          C 0_f1_C

# histogram 
hist(df$richness, main = "histogram of richness")
(f <- ggplot(df, aes(x = Day, y = richness, shape = Frequency, color = Soil)) + 
    geom_point(size = 2, alpha=0.9) +
    labs(x = "Days", y = "Richness")+
    theme_light())

(f <- ggplot(df, aes(x = Day, y = richness, shape = Frequency, color = Soil)) + 
    geom_point(size = 2, alpha=0.9) +
    facet_grid(Soil ~ Frequency) +
    labs(x = "Days", y = "Richness")+
    theme_light())

(f <- ggplot(df, aes(x = factor(Day), y = richness, shape = Frequency, color = Soil)) + 
    geom_boxplot() +
    geom_point(size = 2, alpha=0.9, position = position_jitterdodge()) +
    labs(x = "Days", y = "Richness")+
    theme_light())

1.2 exploring analysis

1.2.1 Generalized linear model

• Richness is a discrete varaible and belongs to discrete distributions, such as poisson, binomial, negbin. Normal, amma, beta distributions are continuous distribution.
• Significant test and q-q plot showed that richness is not or marginally normal distributed.

# full model
glm0 <- glm(richness~Soil*Day*Frequency, family="poisson", data = df)
summary(glm0)

## 
## Call:
## glm(formula = richness ~ Soil * Day * Frequency, family = "poisson", 
##     data = df)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -17.3265   -1.4315    0.2237    2.4846    8.3313  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             6.686904   0.015134 441.859  < 2e-16 ***
## Soil70                  0.340859   0.019967  17.071  < 2e-16 ***
## Day                     0.017359   0.001186  14.637  < 2e-16 ***
## Frequencyf2             0.163617   0.020703   7.903 2.72e-15 ***
## Frequencyf3             0.131680   0.020845   6.317 2.67e-10 ***
## Frequencyf4             0.068546   0.021003   3.264 0.001100 ** 
## Soil70:Day             -0.007455   0.001579  -4.721 2.34e-06 ***
## Soil70:Frequencyf2     -0.124594   0.027741  -4.491 7.08e-06 ***
## Soil70:Frequencyf3     -0.104908   0.027846  -3.767 0.000165 ***
## Soil70:Frequencyf4      0.122734   0.027641   4.440 8.99e-06 ***
## Day:Frequencyf2        -0.004736   0.001633  -2.901 0.003722 ** 
## Day:Frequencyf3        -0.004244   0.001643  -2.583 0.009783 ** 
## Day:Frequencyf4         0.001474   0.001643   0.898 0.369436    
## Soil70:Day:Frequencyf2 -0.004220   0.002215  -1.905 0.056825 .  
## Soil70:Day:Frequencyf3 -0.002340   0.002219  -1.055 0.291584    
## Soil70:Day:Frequencyf4 -0.018084   0.002206  -8.196 2.48e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4746.1  on 143  degrees of freedom
## Residual deviance: 2489.2  on 128  degrees of freedom
## AIC: 3796.2
## 
## Number of Fisher Scoring iterations: 4

plot(glm0)

# check the mean and var of Richness
mean(df$richness); var(df$richness)

## [1] 1131.972

## [1] 36103.89

# check overdispersion
dispersiontest(glm0)

## 
##  Overdispersion test
## 
## data:  glm0
## z = 5.0738, p-value = 1.95e-07
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion 
##   16.18844

# update full model
glm1 <- glm(richness ~ Soil*Day*Frequency, family = "quasipoisson", data = df)
summary(glm1)

## 
## Call:
## glm(formula = richness ~ Soil * Day * Frequency, family = "quasipoisson", 
##     data = df)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -17.3265   -1.4315    0.2237    2.4846    8.3313  
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             6.686904   0.064583 103.539  < 2e-16 ***
## Soil70                  0.340859   0.085211   4.000 0.000106 ***
## Day                     0.017359   0.005061   3.430 0.000813 ***
## Frequencyf2             0.163617   0.088350   1.852 0.066343 .  
## Frequencyf3             0.131680   0.088959   1.480 0.141267    
## Frequencyf4             0.068546   0.089633   0.765 0.445835    
## Soil70:Day             -0.007455   0.006738  -1.106 0.270660    
## Soil70:Frequencyf2     -0.124594   0.118385  -1.052 0.294577    
## Soil70:Frequencyf3     -0.104908   0.118833  -0.883 0.378987    
## Soil70:Frequencyf4      0.122734   0.117961   1.040 0.300087    
## Day:Frequencyf2        -0.004736   0.006967  -0.680 0.497894    
## Day:Frequencyf3        -0.004244   0.007011  -0.605 0.546012    
## Day:Frequencyf4         0.001474   0.007010   0.210 0.833756    
## Soil70:Day:Frequencyf2 -0.004220   0.009454  -0.446 0.656126    
## Soil70:Day:Frequencyf3 -0.002340   0.009471  -0.247 0.805201    
## Soil70:Day:Frequencyf4 -0.018084   0.009416  -1.921 0.057010 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 18.21199)
## 
##     Null deviance: 4746.1  on 143  degrees of freedom
## Residual deviance: 2489.2  on 128  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

# update model by removing non-significant factors
glm2 <- glm(richness ~ Soil * Day, family = "quasipoisson", data = df)
summary(glm2)

## 
## Call:
## glm(formula = richness ~ Soil * Day, family = "quasipoisson", 
##     data = df)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -18.6252   -1.6266    0.4184    2.8346    7.7032  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.779524   0.031647 214.226  < 2e-16 ***
## Soil70       0.313968   0.042300   7.422 1.02e-11 ***
## Day          0.015424   0.002493   6.186 6.39e-09 ***
## Soil70:Day  -0.013585   0.003392  -4.005    1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 18.98598)
## 
##     Null deviance: 4746.1  on 143  degrees of freedom
## Residual deviance: 2851.6  on 140  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

1.2.2 applying general linear model to log transformed richness

# full model
lm0 <- lm(log(richness) ~ Soil * Day * Frequency, data=df)
plot(lm0)
plot(lm0$residuals~df$Day)
plot(lm0$residuals~df$Soil)
plot(lm0$residuals~df$Frequency)
summary(lm0)

## 
## Call:
## lm(formula = log(richness) ~ Soil * Day * Frequency, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.58917 -0.03811  0.01901  0.07819  0.23875 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             6.6708410  0.0638063 104.548  < 2e-16 ***
## Soil70                  0.3309207  0.0902357   3.667 0.000358 ***
## Day                     0.0180250  0.0052567   3.429 0.000816 ***
## Frequencyf2             0.1808577  0.0902357   2.004 0.047150 *  
## Frequencyf3             0.1341841  0.0902357   1.487 0.139464    
## Frequencyf4             0.0813477  0.0902357   0.902 0.369015    
## Soil70:Day             -0.0069091  0.0074341  -0.929 0.354441    
## Soil70:Frequencyf2     -0.1202914  0.1276126  -0.943 0.347646    
## Soil70:Frequencyf3     -0.1041991  0.1276126  -0.817 0.415716    
## Soil70:Frequencyf4      0.1322008  0.1276126   1.036 0.302177    
## Day:Frequencyf2        -0.0056138  0.0074341  -0.755 0.451554    
## Day:Frequencyf3        -0.0041081  0.0074341  -0.553 0.581495    
## Day:Frequencyf4         0.0006347  0.0074341   0.085 0.932092    
## Soil70:Day:Frequencyf2 -0.0057465  0.0105134  -0.547 0.585616    
## Soil70:Day:Frequencyf3 -0.0025173  0.0105134  -0.239 0.811148    
## Soil70:Day:Frequencyf4 -0.0185600  0.0105134  -1.765 0.079887 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.142 on 128 degrees of freedom
## Multiple R-squared:  0.4508, Adjusted R-squared:  0.3864 
## F-statistic: 7.004 on 15 and 128 DF,  p-value: 5.976e-11

# update full model
lm0 <- lm(log(richness) ~ Soil * Day, data=df)
plot(lm0)
summary(lm0)

## 
## Call:
## lm(formula = log(richness) ~ Soil * Day, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.64275 -0.03683  0.01994  0.09358  0.22086 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.769938   0.032648 207.360  < 2e-16 ***
## Soil70       0.307848   0.046171   6.668 5.55e-10 ***
## Day          0.015753   0.002690   5.857 3.22e-08 ***
## Soil70:Day  -0.013615   0.003804  -3.579 0.000474 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1453 on 140 degrees of freedom
## Multiple R-squared:  0.3709, Adjusted R-squared:  0.3575 
## F-statistic: 27.52 on 3 and 140 DF,  p-value: 4.742e-14

1.2.3 linear mixed-effects models (suggested by Hou Meng)

# first approach ---
# full model
lmm0_1 <- lme(log(richness) ~ Soil * Day *Frequency, random=~1|box, data=df, correlation=corAR1(form = ~Day|box))
summary(lmm0_1)

## Linear mixed-effects model fit by REML
##  Data: df 
##        AIC      BIC   logLik
##   -22.5273 31.66128 30.26365
## 
## Random effects:
##  Formula: ~1 | box
##          (Intercept)  Residual
## StdDev: 3.786054e-06 0.1420281
## 
## Correlation Structure: ARMA(1,0)
##  Formula: ~Day | box 
##  Parameter estimate(s):
## Phi1 
##    0 
## Fixed effects: log(richness) ~ Soil * Day * Frequency 
##                            Value  Std.Error  DF   t-value p-value
## (Intercept)             6.670841 0.06380629 112 104.54833  0.0000
## Soil70                  0.330921 0.09023572  16   3.66729  0.0021
## Day                     0.018025 0.00525670 112   3.42895  0.0008
## Frequencyf2             0.180858 0.09023572  16   2.00428  0.0623
## Frequencyf3             0.134184 0.09023572  16   1.48704  0.1564
## Frequencyf4             0.081348 0.09023572  16   0.90150  0.3807
## Soil70:Day             -0.006909 0.00743409 112  -0.92938  0.3547
## Soil70:Frequencyf2     -0.120291 0.12761258  16  -0.94263  0.3599
## Soil70:Frequencyf3     -0.104199 0.12761258  16  -0.81653  0.4262
## Soil70:Frequencyf4      0.132201 0.12761258  16   1.03595  0.3156
## Day:Frequencyf2        -0.005614 0.00743409 112  -0.75514  0.4518
## Day:Frequencyf3        -0.004108 0.00743409 112  -0.55261  0.5816
## Day:Frequencyf4         0.000635 0.00743409 112   0.08538  0.9321
## Soil70:Day:Frequencyf2 -0.005746 0.01051339 112  -0.54658  0.5858
## Soil70:Day:Frequencyf3 -0.002517 0.01051339 112  -0.23944  0.8112
## Soil70:Day:Frequencyf4 -0.018560 0.01051339 112  -1.76536  0.0802
##  Correlation: 
##                        (Intr) Soil70 Day    Frqnc2 Frqnc3 Frqnc4 Sl70:D S70:F2
## Soil70                 -0.707                                                 
## Day                    -0.851  0.602                                          
## Frequencyf2            -0.707  0.500  0.602                                   
## Frequencyf3            -0.707  0.500  0.602  0.500                            
## Frequencyf4            -0.707  0.500  0.602  0.500  0.500                     
## Soil70:Day              0.602 -0.851 -0.707 -0.426 -0.426 -0.426              
## Soil70:Frequencyf2      0.500 -0.707 -0.426 -0.707 -0.354 -0.354  0.602       
## Soil70:Frequencyf3      0.500 -0.707 -0.426 -0.354 -0.707 -0.354  0.602  0.500
## Soil70:Frequencyf4      0.500 -0.707 -0.426 -0.354 -0.354 -0.707  0.602  0.500
## Day:Frequencyf2         0.602 -0.426 -0.707 -0.851 -0.426 -0.426  0.500  0.602
## Day:Frequencyf3         0.602 -0.426 -0.707 -0.426 -0.851 -0.426  0.500  0.301
## Day:Frequencyf4         0.602 -0.426 -0.707 -0.426 -0.426 -0.851  0.500  0.301
## Soil70:Day:Frequencyf2 -0.426  0.602  0.500  0.602  0.301  0.301 -0.707 -0.851
## Soil70:Day:Frequencyf3 -0.426  0.602  0.500  0.301  0.602  0.301 -0.707 -0.426
## Soil70:Day:Frequencyf4 -0.426  0.602  0.500  0.301  0.301  0.602 -0.707 -0.426
##                        S70:F3 S70:F4 Dy:Fr2 Dy:Fr3 Dy:Fr4 S70:D:F2 S70:D:F3
## Soil70                                                                     
## Day                                                                        
## Frequencyf2                                                                
## Frequencyf3                                                                
## Frequencyf4                                                                
## Soil70:Day                                                                 
## Soil70:Frequencyf2                                                         
## Soil70:Frequencyf3                                                         
## Soil70:Frequencyf4      0.500                                              
## Day:Frequencyf2         0.301  0.301                                       
## Day:Frequencyf3         0.602  0.301  0.500                                
## Day:Frequencyf4         0.301  0.602  0.500  0.500                         
## Soil70:Day:Frequencyf2 -0.426 -0.426 -0.707 -0.354 -0.354                  
## Soil70:Day:Frequencyf3 -0.851 -0.426 -0.354 -0.707 -0.354  0.500           
## Soil70:Day:Frequencyf4 -0.426 -0.851 -0.354 -0.354 -0.707  0.500    0.500  
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -4.1482469 -0.2683141  0.1338492  0.5505475  1.6809933 
## 
## Number of Observations: 144
## Number of Groups: 24

# update model
lmm1_1 <- lme(log(richness) ~ Soil * Day, random=~1|box, data=df, correlation=corAR1(form = ~Day|box))
summary(lmm1_1)

## Linear mixed-effects model fit by REML
##  Data: df 
##         AIC       BIC   logLik
##   -104.2061 -83.61463 59.10306
## 
## Random effects:
##  Formula: ~1 | box
##          (Intercept)  Residual
## StdDev: 1.176309e-05 0.1453448
## 
## Correlation Structure: ARMA(1,0)
##  Formula: ~Day | box 
##  Parameter estimate(s):
## Phi1 
##    0 
## Fixed effects: log(richness) ~ Soil * Day 
##                 Value  Std.Error  DF   t-value p-value
## (Intercept)  6.769938 0.03264816 118 207.36050   0e+00
## Soil70       0.307848 0.04617147  22   6.66750   0e+00
## Day          0.015753 0.00268973 118   5.85679   0e+00
## Soil70:Day  -0.013615 0.00380385 118  -3.57928   5e-04
##  Correlation: 
##            (Intr) Soil70 Day   
## Soil70     -0.707              
## Day        -0.851  0.602       
## Soil70:Day  0.602 -0.851 -0.707
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -4.4222157 -0.2534009  0.1372206  0.6438246  1.5195804 
## 
## Number of Observations: 144
## Number of Groups: 24

# second approach ---
# full model
lmm0_2 <- lme(log(richness) ~ Soil * Day *Frequency, random=~1|box, data=df)
summary(lmm0_2)

## Linear mixed-effects model fit by REML
##  Data: df 
##        AIC      BIC   logLik
##   -24.5273 26.80925 30.26365
## 
## Random effects:
##  Formula: ~1 | box
##          (Intercept)  Residual
## StdDev: 3.785982e-06 0.1420281
## 
## Fixed effects: log(richness) ~ Soil * Day * Frequency 
##                            Value  Std.Error  DF   t-value p-value
## (Intercept)             6.670841 0.06380629 112 104.54833  0.0000
## Soil70                  0.330921 0.09023572  16   3.66729  0.0021
## Day                     0.018025 0.00525670 112   3.42895  0.0008
## Frequencyf2             0.180858 0.09023572  16   2.00428  0.0623
## Frequencyf3             0.134184 0.09023572  16   1.48704  0.1564
## Frequencyf4             0.081348 0.09023572  16   0.90150  0.3807
## Soil70:Day             -0.006909 0.00743409 112  -0.92938  0.3547
## Soil70:Frequencyf2     -0.120291 0.12761258  16  -0.94263  0.3599
## Soil70:Frequencyf3     -0.104199 0.12761258  16  -0.81653  0.4262
## Soil70:Frequencyf4      0.132201 0.12761258  16   1.03595  0.3156
## Day:Frequencyf2        -0.005614 0.00743409 112  -0.75514  0.4518
## Day:Frequencyf3        -0.004108 0.00743409 112  -0.55261  0.5816
## Day:Frequencyf4         0.000635 0.00743409 112   0.08538  0.9321
## Soil70:Day:Frequencyf2 -0.005746 0.01051339 112  -0.54658  0.5858
## Soil70:Day:Frequencyf3 -0.002517 0.01051339 112  -0.23944  0.8112
## Soil70:Day:Frequencyf4 -0.018560 0.01051339 112  -1.76536  0.0802
##  Correlation: 
##                        (Intr) Soil70 Day    Frqnc2 Frqnc3 Frqnc4 Sl70:D S70:F2
## Soil70                 -0.707                                                 
## Day                    -0.851  0.602                                          
## Frequencyf2            -0.707  0.500  0.602                                   
## Frequencyf3            -0.707  0.500  0.602  0.500                            
## Frequencyf4            -0.707  0.500  0.602  0.500  0.500                     
## Soil70:Day              0.602 -0.851 -0.707 -0.426 -0.426 -0.426              
## Soil70:Frequencyf2      0.500 -0.707 -0.426 -0.707 -0.354 -0.354  0.602       
## Soil70:Frequencyf3      0.500 -0.707 -0.426 -0.354 -0.707 -0.354  0.602  0.500
## Soil70:Frequencyf4      0.500 -0.707 -0.426 -0.354 -0.354 -0.707  0.602  0.500
## Day:Frequencyf2         0.602 -0.426 -0.707 -0.851 -0.426 -0.426  0.500  0.602
## Day:Frequencyf3         0.602 -0.426 -0.707 -0.426 -0.851 -0.426  0.500  0.301
## Day:Frequencyf4         0.602 -0.426 -0.707 -0.426 -0.426 -0.851  0.500  0.301
## Soil70:Day:Frequencyf2 -0.426  0.602  0.500  0.602  0.301  0.301 -0.707 -0.851
## Soil70:Day:Frequencyf3 -0.426  0.602  0.500  0.301  0.602  0.301 -0.707 -0.426
## Soil70:Day:Frequencyf4 -0.426  0.602  0.500  0.301  0.301  0.602 -0.707 -0.426
##                        S70:F3 S70:F4 Dy:Fr2 Dy:Fr3 Dy:Fr4 S70:D:F2 S70:D:F3
## Soil70                                                                     
## Day                                                                        
## Frequencyf2                                                                
## Frequencyf3                                                                
## Frequencyf4                                                                
## Soil70:Day                                                                 
## Soil70:Frequencyf2                                                         
## Soil70:Frequencyf3                                                         
## Soil70:Frequencyf4      0.500                                              
## Day:Frequencyf2         0.301  0.301                                       
## Day:Frequencyf3         0.602  0.301  0.500                                
## Day:Frequencyf4         0.301  0.602  0.500  0.500                         
## Soil70:Day:Frequencyf2 -0.426 -0.426 -0.707 -0.354 -0.354                  
## Soil70:Day:Frequencyf3 -0.851 -0.426 -0.354 -0.707 -0.354  0.500           
## Soil70:Day:Frequencyf4 -0.426 -0.851 -0.354 -0.354 -0.707  0.500    0.500  
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -4.1482469 -0.2683141  0.1338492  0.5505475  1.6809933 
## 
## Number of Observations: 144
## Number of Groups: 24

# update model
lmm1_2 <- lme(log(richness) ~ Soil * Day, random=~1|box, data=df)
summary(lmm1_2)

## Linear mixed-effects model fit by REML
##  Data: df 
##         AIC       BIC   logLik
##   -106.2061 -88.55627 59.10306
## 
## Random effects:
##  Formula: ~1 | box
##          (Intercept)  Residual
## StdDev: 1.134264e-05 0.1453448
## 
## Fixed effects: log(richness) ~ Soil * Day 
##                 Value  Std.Error  DF   t-value p-value
## (Intercept)  6.769938 0.03264816 118 207.36050   0e+00
## Soil70       0.307848 0.04617147  22   6.66750   0e+00
## Day          0.015753 0.00268973 118   5.85679   0e+00
## Soil70:Day  -0.013615 0.00380385 118  -3.57928   5e-04
##  Correlation: 
##            (Intr) Soil70 Day   
## Soil70     -0.707              
## Day        -0.851  0.602       
## Soil70:Day  0.602 -0.851 -0.707
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -4.4222157 -0.2534009  0.1372206  0.6438246  1.5195804 
## 
## Number of Observations: 144
## Number of Groups: 24

1.2.4 Generalized least squares fit by REML

# full model which considers non homogeneouse of variance
gls0 <- gls(log(richness) ~ Soil * Day * Frequency, correlation=corAR1(form= ~Day |box), weights=varIdent(form = ~ 1|Soil), na.action=na.omit, data=df)
summary(gls0)

## Generalized least squares fit by REML
##   Model: log(richness) ~ Soil * Day * Frequency 
##   Data: df 
##         AIC      BIC   logLik
##   -33.79917 20.38941 35.89958
## 
## Correlation Structure: ARMA(1,0)
##  Formula: ~Day | box 
##  Parameter estimate(s):
## Phi1 
##    0 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Soil 
##  Parameter estimates:
##        0       70 
## 1.000000 1.530899 
## 
## Coefficients:
##                            Value  Std.Error   t-value p-value
## (Intercept)             6.670841 0.04934780 135.18010  0.0000
## Soil70                  0.330921 0.09023571   3.66729  0.0004
## Day                     0.018025 0.00406553   4.43360  0.0000
## Frequencyf2             0.180858 0.06978833   2.59152  0.0107
## Frequencyf3             0.134184 0.06978833   1.92273  0.0567
## Frequencyf4             0.081348 0.06978833   1.16563  0.2459
## Soil70:Day             -0.006909 0.00743409  -0.92938  0.3544
## Soil70:Frequencyf2     -0.120291 0.12761257  -0.94263  0.3476
## Soil70:Frequencyf3     -0.104199 0.12761257  -0.81653  0.4157
## Soil70:Frequencyf4      0.132201 0.12761257   1.03595  0.3022
## Day:Frequencyf2        -0.005614 0.00574953  -0.97639  0.3307
## Day:Frequencyf3        -0.004108 0.00574953  -0.71452  0.4762
## Day:Frequencyf4         0.000635 0.00574953   0.11040  0.9123
## Soil70:Day:Frequencyf2 -0.005746 0.01051339  -0.54658  0.5856
## Soil70:Day:Frequencyf3 -0.002517 0.01051339  -0.23944  0.8111
## Soil70:Day:Frequencyf4 -0.018560 0.01051339  -1.76536  0.0799
## 
##  Correlation: 
##                        (Intr) Soil70 Day    Frqnc2 Frqnc3 Frqnc4 Sl70:D S70:F2
## Soil70                 -0.547                                                 
## Day                    -0.851  0.466                                          
## Frequencyf2            -0.707  0.387  0.602                                   
## Frequencyf3            -0.707  0.387  0.602  0.500                            
## Frequencyf4            -0.707  0.387  0.602  0.500  0.500                     
## Soil70:Day              0.466 -0.851 -0.547 -0.329 -0.329 -0.329              
## Soil70:Frequencyf2      0.387 -0.707 -0.329 -0.547 -0.273 -0.273  0.602       
## Soil70:Frequencyf3      0.387 -0.707 -0.329 -0.273 -0.547 -0.273  0.602  0.500
## Soil70:Frequencyf4      0.387 -0.707 -0.329 -0.273 -0.273 -0.547  0.602  0.500
## Day:Frequencyf2         0.602 -0.329 -0.707 -0.851 -0.426 -0.426  0.387  0.466
## Day:Frequencyf3         0.602 -0.329 -0.707 -0.426 -0.851 -0.426  0.387  0.233
## Day:Frequencyf4         0.602 -0.329 -0.707 -0.426 -0.426 -0.851  0.387  0.233
## Soil70:Day:Frequencyf2 -0.329  0.602  0.387  0.466  0.233  0.233 -0.707 -0.851
## Soil70:Day:Frequencyf3 -0.329  0.602  0.387  0.233  0.466  0.233 -0.707 -0.426
## Soil70:Day:Frequencyf4 -0.329  0.602  0.387  0.233  0.233  0.466 -0.707 -0.426
##                        S70:F3 S70:F4 Dy:Fr2 Dy:Fr3 Dy:Fr4 S70:D:F2 S70:D:F3
## Soil70                                                                     
## Day                                                                        
## Frequencyf2                                                                
## Frequencyf3                                                                
## Frequencyf4                                                                
## Soil70:Day                                                                 
## Soil70:Frequencyf2                                                         
## Soil70:Frequencyf3                                                         
## Soil70:Frequencyf4      0.500                                              
## Day:Frequencyf2         0.233  0.233                                       
## Day:Frequencyf3         0.466  0.233  0.500                                
## Day:Frequencyf4         0.233  0.466  0.500  0.500                         
## Soil70:Day:Frequencyf2 -0.426 -0.426 -0.547 -0.273 -0.273                  
## Soil70:Day:Frequencyf3 -0.851 -0.426 -0.273 -0.547 -0.273  0.500           
## Soil70:Day:Frequencyf4 -0.426 -0.851 -0.273 -0.273 -0.547  0.500    0.500  
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.6858047 -0.3073183  0.1130485  0.5706062  1.7771310 
## 
## Residual standard error: 0.1098446 
## Degrees of freedom: 144 total; 128 residual

1.2.5 glmm

# full model
glmm1 <- glmer(richness ~ Soil * Day * Frequency + (Day | box), family = poisson, df)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0115456 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

summary(glmm1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: richness ~ Soil * Day * Frequency + (Day | box)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##   3546.2   3602.6  -1754.1   3508.2      125 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -14.2555  -1.2237   0.2777   2.1606   9.9149 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev. Corr 
##  box    (Intercept) 5.994e-03 0.077424      
##         Day         3.195e-05 0.005652 -0.93
## Number of obs: 144, groups:  box, 24
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             6.681559   0.047209 141.531  < 2e-16 ***
## Soil70                  0.345403   0.066308   5.209 1.90e-07 ***
## Day                     0.017647   0.003473   5.081 3.75e-07 ***
## Frequencyf2             0.168198   0.066532   2.528  0.01147 *  
## Frequencyf3             0.135751   0.066578   2.039  0.04145 *  
## Frequencyf4             0.072005   0.066628   1.081  0.27983    
## Soil70:Day             -0.007814   0.004878  -1.602  0.10923    
## Soil70:Frequencyf2     -0.129669   0.093618  -1.385  0.16603    
## Soil70:Frequencyf3     -0.110441   0.093651  -1.179  0.23829    
## Soil70:Frequencyf4      0.119978   0.093589   1.282  0.19985    
## Day:Frequencyf2        -0.004996   0.004896  -1.020  0.30751    
## Day:Frequencyf3        -0.004522   0.004899  -0.923  0.35602    
## Day:Frequencyf4         0.001281   0.004899   0.262  0.79371    
## Soil70:Day:Frequencyf2 -0.003891   0.006893  -0.565  0.57239    
## Soil70:Day:Frequencyf3 -0.001878   0.006894  -0.272  0.78532    
## Soil70:Day:Frequencyf4 -0.017851   0.006890  -2.591  0.00958 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

## convergence code: 0
## Model failed to converge with max|grad| = 0.0115456 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

# update model
# update model by considering overdispersion
# add a observation-level random effect (sample) if overdispersed Elston et al. (2001).
df$sample <- factor(paste(df$Soil, df$Frequency, df$Day, df$replicates, sep = "_"))
str(df)

## 'data.frame':    144 obs. of  7 variables:
##  $ richness  : int  828 1075 1101 1100 1108 1124 792 900 727 1074 ...
##  $ Soil      : Factor w/ 2 levels "0","70": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Frequency : Factor w/ 4 levels "f1","f2","f3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Day       : int  12 12 12 16 16 16 2 2 2 20 ...
##  $ replicates: Factor w/ 3 levels "A","B","C": 1 2 3 1 2 3 1 2 3 1 ...
##  $ box       : Factor w/ 24 levels "0_f1_A","0_f1_B",..: 1 2 3 1 2 3 1 2 3 1 ...
##  $ sample    : Factor w/ 144 levels "0_f1_12_A","0_f1_12_B",..: 1 2 3 4 5 6 10 11 12 7 ...

glmm2 <- glmer(richness ~ Soil * Day + (Day | box) + (1 | sample), family = poisson, df)

## boundary (singular) fit: see ?isSingular

# check overdispersion https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#testing-for-overdispersioncomputing-overdispersion-factor
overdisp_fun <- function(model) {
  rdf <- df.residual(model)
  rp <- residuals(model,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
# p < 0.05, overdispersed
overdisp_fun(glmm2) 

##        chisq        ratio          rdf            p 
##   8.46045473   0.06220923 136.00000000   1.00000000

tab_model(glmm2,show.icc=T)

## boundary (singular) fit: see ?isSingular

	richness
Predictors	Incidence Rate Ratios	CI	p
(Intercept)	871.77	815.12 – 932.34	<0.001
Soil [70]	1.36	1.24 – 1.50	<0.001
Day	1.02	1.01 – 1.02	<0.001
Soil [70] * Day	0.99	0.98 – 0.99	<0.001
Random Effects
σ2	0.02
τ00 sample	0.02
τ00 box	0.00
τ11 box.Day	0.00
ρ01 box	-1.00
N box	24
N sample	144
Observations	144
Marginal R2 / Conditional R2	0.389 / NA

summary(glmm2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: richness ~ Soil * Day + (Day | box) + (1 | sample)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##   1883.5   1907.3   -933.8   1867.5      136 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.24534 -0.05996  0.02397  0.12886  0.28423 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev. Corr 
##  sample (Intercept) 1.837e-02 0.135532      
##  box    (Intercept) 2.333e-03 0.048301      
##         Day         7.295e-06 0.002701 -1.00
## Number of obs: 144, groups:  sample, 144; box, 24
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  6.770520   0.034277 197.525  < 2e-16 ***
## Soil70       0.308498   0.048363   6.379 1.78e-10 ***
## Day          0.015717   0.002689   5.845 5.07e-09 ***
## Soil70:Day  -0.013607   0.003797  -3.584 0.000339 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) Soil70 Day   
## Soil70     -0.709              
## Day        -0.863  0.611       
## Soil70:Day  0.611 -0.863 -0.708
## convergence code: 0
## boundary (singular) fit: see ?isSingular

anova(glmm2)

## Analysis of Variance Table
##          npar Sum Sq Mean Sq F value
## Soil        1 42.122  42.122  42.122
## Day         1 21.938  21.938  21.938
## Soil:Day    1 12.844  12.844  12.844

glmm的syntax写的是否合适，能否体现时间自相关？ summary结果里的Correlation of Fixed Effects很费解啊。。。

1.2.6 swtich from frequency based appraoch to bayesian based approach (do not work!)

# brm0 <- brm(bf(richness ~ Soil + Day + (Day|box), sigma ~ Soil), data = df, autocor = cor_cosy(~Day|box))
# brm1 <- brm(richness ~ Soil * Day, data = df, family = negbinomial("log"))
# summary(brm1)

1.3 final model - Fit Linear Model to Richness using Generalized Least Squares

1.3.1 full model

# update model by removing not significant factors
gls1 <- gls(log(richness) ~ Soil + Day + Soil:Day, correlation=corAR1(form= ~Day |box), weights=varIdent(form = ~ 1|Soil), na.action=na.omit, data=df)
summary(gls1)

## Generalized least squares fit by REML
##   Model: log(richness) ~ Soil + Day + Soil:Day 
##   Data: df 
##         AIC       BIC   logLik
##   -114.0469 -93.45544 64.02347
## 
## Correlation Structure: ARMA(1,0)
##  Formula: ~Day | box 
##  Parameter estimate(s):
## Phi1 
##    0 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Soil 
##  Parameter estimates:
##        0       70 
## 1.000000 1.461337 
## 
## Coefficients:
##                 Value  Std.Error   t-value p-value
## (Intercept)  6.769938 0.02607474 259.63593   0e+00
## Soil70       0.307848 0.04617147   6.66750   0e+00
## Day          0.015753 0.00214817   7.33328   0e+00
## Soil70:Day  -0.013615 0.00380385  -3.57928   5e-04
## 
##  Correlation: 
##            (Intr) Soil70 Day   
## Soil70     -0.565              
## Day        -0.851  0.481       
## Soil70:Day  0.481 -0.851 -0.565
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.1849151 -0.2226717  0.1280295  0.6379215  1.6976512 
## 
## Residual standard error: 0.1160809 
## Degrees of freedom: 144 total; 140 residual

anova(gls1)

## Denom. DF: 140 
##             numDF  F-value p-value
## (Intercept)     1 382894.3  <.0001
## Soil            1     47.6  <.0001
## Day             1     41.4  <.0001
## Soil:Day        1     12.8   5e-04

plot(gls1)

on day 0, the richness of 0 year soil inundated is \(e^{6.764}\), while the richness of 70 year soil at f1 is \(e^{(6.764+0.331)}\), i.e. \(e^{7.095}\). For the richness of 0 and 70 year soils are:

\(richness_{0yr} = e^{(6.764 + 0.0163*Day)}\); and

\(richness_{70yr} = e^{(7.095 + (0.0163-0.0146)*Day)}\), i.e. \(richness_{70yr} = e^{(7.095 + 0.0018*Day)}\), respectively.

1.3.2 model diagnosis

1.1.3.1 compare with initial model

# compare two models
AIC(gls0, gls1)

## Warning in AIC.default(gls0, gls1): models are not all fitted to the same number
## of observations

##      df        AIC
## gls0 19  -33.79917
## gls1  7 -114.04694

plot(gls1)

# check normality via histogram (1)
hist(gls1$residuals, main="Residual distribution")

## assess the normality of residuals using shapiro wilk test (2)
shapiro.test(gls1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  gls1$residuals
## W = 0.84583, p-value = 5.526e-11

plot(gls1$residuals ~ gls1$fitted)
(f <- ggplot(df, aes(x = Day, y = log(richness), shape = Frequency, color = Soil)) + 
    geom_point(size = 2, alpha=0.9) +
    facet_grid(Soil ~ Frequency) +
    labs(x = "Days", y = "Richness")+
    theme_light())

the Shapiro Wilk test was significant (p < 0.05), indicating residuals are not normal distributed?!

2. both sea water and strile water treatments

2.1 load dataset (2 soil types x 2 water treatments x 2 inundation frequency)

# read dataset
df <- read.csv("alpha_diversity_dispersal.csv", row.names = 1, check.names = 0)
# only keep richness
df <- df[, -c(2:7)]

# remove samples from day 0
df <- subset(df, !grepl('___', rownames(df)))

# select soil from sea water treatment
df <- df[df$Frequency == "f1" | df$Frequency == "f4",]
df$Soil <- factor(df$Soil, levels = c("0", "70"))
df$Water <- factor(df$Water)
df$Frequency <- factor(df$Frequency)
df$box <- factor(paste(df$Soil, df$Frequency, df$replicates, sep = "_"))
str(df)

## 'data.frame':    144 obs. of  7 variables:
##  $ richness  : int  828 1075 1101 1100 1108 1124 792 900 727 1074 ...
##  $ Soil      : Factor w/ 2 levels "0","70": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Frequency : Factor w/ 2 levels "f1","f4": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Water     : Factor w/ 2 levels "sea","strl": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Day       : int  12 12 12 16 16 16 2 2 2 20 ...
##  $ replicates: Factor w/ 3 levels "A","B","C": 1 2 3 1 2 3 1 2 3 1 ...
##  $ box       : Factor w/ 12 levels "0_f1_A","0_f1_B",..: 1 2 3 1 2 3 1 2 3 1 ...

head(df)

##               richness Soil Frequency Water Day replicates    box
## 0_f1_sea_12_A      828    0        f1   sea  12          A 0_f1_A
## 0_f1_sea_12_B     1075    0        f1   sea  12          B 0_f1_B
## 0_f1_sea_12_C     1101    0        f1   sea  12          C 0_f1_C
## 0_f1_sea_16_A     1100    0        f1   sea  16          A 0_f1_A
## 0_f1_sea_16_B     1108    0        f1   sea  16          B 0_f1_B
## 0_f1_sea_16_C     1124    0        f1   sea  16          C 0_f1_C

# histogram 
hist(df$richness, main = "histogram of richness")
(f <- ggplot(df, aes(x = Day, y = richness, shape = Frequency, color = Water)) + 
    geom_point(size = 2, alpha=0.9) +
    facet_grid(Soil ~ Frequency) +
    labs(x = "Days", y = "Richness")+
    theme_light())

(f <- ggplot(df, aes(x = factor(Day), y = richness, shape = Frequency, color = Water)) +
    geom_boxplot() +
    facet_grid(Soil ~ .) +
    geom_point(size = 2, alpha=0.9, position = position_jitterdodge()) +
    labs(x = "Days", y = "Richness")+
    theme_light())

2.2 Fit Linear Model to Richness using Generalized Least Squares

2.2.1 full model

# full model 

gls0 <- gls(log(richness) ~ Soil * Day * Water * Frequency, correlation=corAR1(form= ~1 |box), weights=varIdent(form = ~ 1|Soil), na.action=na.omit, data=df)
summary(gls0)

## Generalized least squares fit by REML
##   Model: log(richness) ~ Soil * Day * Water * Frequency 
##   Data: df 
##         AIC      BIC   logLik
##   -34.53114 19.65744 36.26557
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | box 
##  Parameter estimate(s):
##         Phi 
## -0.04157485 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Soil 
##  Parameter estimates:
##         0        70 
## 1.0000000 0.7923043 
## 
## Coefficients:
##                                      Value  Std.Error  t-value p-value
## (Intercept)                       6.670819 0.06919295 96.40895  0.0000
## Soil70                            0.330203 0.08827854  3.74047  0.0003
## Day                               0.018140 0.00579290  3.13139  0.0022
## Waterstrl                         0.150218 0.09792631  1.53399  0.1275
## Frequencyf4                       0.082946 0.09785360  0.84765  0.3982
## Soil70:Day                       -0.007070 0.00739076 -0.95662  0.3406
## Soil70:Waterstrl                 -0.078749 0.12493747 -0.63031  0.5296
## Day:Waterstrl                    -0.011984 0.00818955 -1.46329  0.1458
## Soil70:Frequencyf4                0.131631 0.12484471  1.05436  0.2937
## Day:Frequencyf4                   0.000335 0.00819240  0.04090  0.9674
## Waterstrl:Frequencyf4            -0.215139 0.13848871 -1.55347  0.1228
## Soil70:Day:Waterstrl              0.008981 0.01044849  0.85958  0.3916
## Soil70:Day:Frequencyf4           -0.018248 0.01045212 -1.74583  0.0832
## Soil70:Waterstrl:Frequencyf4      0.032896 0.17668826  0.18618  0.8526
## Day:Waterstrl:Frequencyf4         0.012234 0.01158178  1.05628  0.2928
## Soil70:Day:Waterstrl:Frequencyf4 -0.002956 0.01477640 -0.20003  0.8418
## 
##  Correlation: 
##                                  (Intr) Soil70 Day    Wtrstr Frqnc4 Sl70:D
## Soil70                           -0.784                                   
## Day                              -0.865  0.678                            
## Waterstrl                        -0.708  0.555  0.612                     
## Frequencyf4                      -0.707  0.554  0.612  0.500              
## Soil70:Day                        0.678 -0.865 -0.784 -0.480 -0.480       
## Soil70:Waterstrl                  0.555 -0.708 -0.480 -0.784 -0.392  0.612
## Day:Waterstrl                     0.611 -0.479 -0.707 -0.865 -0.432  0.554
## Soil70:Frequencyf4                0.554 -0.707 -0.480 -0.392 -0.784  0.612
## Day:Frequencyf4                   0.612 -0.480 -0.707 -0.433 -0.865  0.554
## Waterstrl:Frequencyf4             0.500 -0.392 -0.433 -0.707 -0.708  0.339
## Soil70:Day:Waterstrl             -0.479  0.611  0.554  0.678  0.339 -0.707
## Soil70:Day:Frequencyf4           -0.480  0.612  0.554  0.339  0.678 -0.707
## Soil70:Waterstrl:Frequencyf4     -0.392  0.500  0.339  0.554  0.555 -0.433
## Day:Waterstrl:Frequencyf4        -0.432  0.339  0.500  0.611  0.611 -0.392
## Soil70:Day:Waterstrl:Frequencyf4  0.339 -0.432 -0.392 -0.479 -0.479  0.500
##                                  Sl70:W Dy:Wtr S70:F4 Dy:Fr4 Wtr:F4 Sl70:D:W
## Soil70                                                                      
## Day                                                                         
## Waterstrl                                                                   
## Frequencyf4                                                                 
## Soil70:Day                                                                  
## Soil70:Waterstrl                                                            
## Day:Waterstrl                     0.678                                     
## Soil70:Frequencyf4                0.500  0.339                              
## Day:Frequencyf4                   0.339  0.500  0.678                       
## Waterstrl:Frequencyf4             0.554  0.611  0.555  0.612                
## Soil70:Day:Waterstrl             -0.865 -0.784 -0.432 -0.392 -0.479         
## Soil70:Day:Frequencyf4           -0.433 -0.392 -0.865 -0.784 -0.480  0.500  
## Soil70:Waterstrl:Frequencyf4     -0.707 -0.479 -0.708 -0.480 -0.784  0.611  
## Day:Waterstrl:Frequencyf4        -0.479 -0.707 -0.479 -0.707 -0.865  0.554  
## Soil70:Day:Waterstrl:Frequencyf4  0.611  0.554  0.611  0.554  0.678 -0.707  
##                                  S70:D:F S70:W: D:W:F4
## Soil70                                                
## Day                                                   
## Waterstrl                                             
## Frequencyf4                                           
## Soil70:Day                                            
## Soil70:Waterstrl                                      
## Day:Waterstrl                                         
## Soil70:Frequencyf4                                    
## Day:Frequencyf4                                       
## Waterstrl:Frequencyf4                                 
## Soil70:Day:Waterstrl                                  
## Soil70:Day:Frequencyf4                                
## Soil70:Waterstrl:Frequencyf4      0.612               
## Day:Waterstrl:Frequencyf4         0.554   0.678       
## Soil70:Day:Waterstrl:Frequencyf4 -0.707  -0.865 -0.784
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.6712668 -0.4130164  0.1313822  0.6485909  1.4405870 
## 
## Residual standard error: 0.152314 
## Degrees of freedom: 144 total; 128 residual

# update model by removing not significant factors
gls1 <- gls(log(richness) ~ Soil + Day + Soil:Day, correlation=corAR1(form= ~1 |box), weights=varIdent(form = ~ 1|Soil),  na.action=na.omit, data=df)
summary(gls1)

## Generalized least squares fit by REML
##   Model: log(richness) ~ Soil + Day + Soil:Day 
##   Data: df 
##         AIC       BIC   logLik
##   -117.3097 -96.71821 65.65486
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | box 
##  Parameter estimate(s):
##         Phi 
## -0.07361147 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Soil 
##  Parameter estimates:
##         0        70 
## 1.0000000 0.8203663 
## 
## Coefficients:
##                 Value  Std.Error   t-value p-value
## (Intercept)  6.733888 0.03512164 191.73046  0.0000
## Soil70       0.364713 0.04542792   8.02838  0.0000
## Day          0.015352 0.00297846   5.15439  0.0000
## Soil70:Day  -0.012466 0.00385247  -3.23582  0.0015
## 
##  Correlation: 
##            (Intr) Soil70 Day   
## Soil70     -0.773              
## Day        -0.877  0.678       
## Soil70:Day  0.678 -0.877 -0.773
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.9921583 -0.4655684  0.1743520  0.6703851  1.6376582 
## 
## Residual standard error: 0.1534817 
## Degrees of freedom: 144 total; 140 residual

anova(gls1)

## Denom. DF: 140 
##             numDF  F-value p-value
## (Intercept)     1 430638.5  <.0001
## Soil            1    116.4  <.0001
## Day             1     17.5  0.0001
## Soil:Day        1     10.5  0.0015

on day 0, the richness of 0 year soil inundated is \(e^{6.733}\), while the richness of 70 year soil at f1 is \(e^{(6.733+0.365)}\), i.e. \(e^{7.098}\). For the richness of 0 and 70 year soils are:

\(richness_{0yr} = e^{(6.733 + 0.0154*Day)}\); and

\(richness_{70yr} = e^{(7.098 + (0.0154-0.0124)*Day)}\), i.e. \(richness_{70yr} = e^{(7.098 + 0.0030*Day)}\), respectively.

1.3.2 model diagnosis

1.1.3.1 compare with initial model

# compare two models
AIC(gls0, gls1)

## Warning in AIC.default(gls0, gls1): models are not all fitted to the same number
## of observations

##      df        AIC
## gls0 19  -34.53114
## gls1  7 -117.30971

plot(gls1)

# check normality via histogram (1)
hist(gls1$residuals, main="Residual distribution")

## assess the normality of residuals using shapiro wilk test (2)
shapiro.test(gls1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  gls1$residuals
## W = 0.90314, p-value = 3.324e-08

plot(gls1$residuals ~ gls1$fitted)

(f <- ggplot(df, aes(x = Day, y = log(richness), shape = Frequency, color = Soil)) + 
    geom_point(size = 2, alpha=0.9) +
    facet_grid(Soil ~ Frequency) +
    labs(x = "Days", y = "Richness")+
    theme_light())

the Shapiro Wilk test was significant (p < 0.05), indicating residuals are not normal distributed?!